% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

%z_stochss

%% Generate random shocks

options = optimset;
T = 10000;               % simulation periods
t_drop = 1000;

rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock
%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
pol_flag   = zeros(T+1,1);

rhsee_derived_temp = zeros(T+1,1);
rbar_derived       = zeros(T+1,1);
rbar_gap           = zeros(T+1,1);

%% Simulation

pol1_t = 0;
pol2_t = 0;

kv(1)  = k_ss;
rbv(1) = rb_ss;
nv(1)  = n_ss;
lv(1)  = kv(1)/nv(1);
av(1)  = a_ss;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t-1);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x*solved_eta1_rhsee);
                rbv(t)    = exp(x*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x*solved_eta3_rhsee);
                rbv(t)    = exp(x*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x*solved_eta2_rhsee);
            rbv(t)    = exp(x*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk_next = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        % check policy
        if pibv(t) >= thre_pib && lv(t)>thre_l
            
            if nv(t) < thre_n
                lv(t)  = lv(t);
                rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk_next,lambda,gamma,delta,sigma_rk);   % deposit interest rate                    
            elseif nv(t) < (1+pn)*thre_n
                pol1_t = pol1_t + 1;
                lv(t) = lv(t) - (nv(t)-thre_n)/((1+pn)*thre_n-thre_n)*(lv(t)-thre_l);
                kv(t) = lv(t)*nv(t);
                iv(t) = kv(t) - (1-delta)*kv(t-1);
                cv(t) = yv(t) - iv(t);
                rhseev(t) = exp( sum(x*solved_eta4_rhsee ) );
                rbv(t)    = 1/( cv(t) - psi*hv(t)^(1+1/nu)/(1+1/nu) ) / rhseev(t);
            else
                pol2_t = pol2_t + 1;
                lv(t) = thre_l;
                kv(t) = lv(t)*nv(t);
                iv(t) = kv(t) - (1-delta)*kv(t-1);
                cv(t) = yv(t) - iv(t);
                rhseev(t) = exp( sum(x*solved_eta4_rhsee ) );
                rbv(t)    = 1/( cv(t) - psi*hv(t)^(1+1/nu)/(1+1/nu) ) / rhseev(t);
            end

        else
            rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk_next,lambda,gamma,delta,sigma_rk);   % deposit interest rate                    
        end

        % run case
        if xiv(t)<1
            lv(t) = thre_l;
            kv(t) = lv(t)*nv(t);
            iv(t) = kv(t) - (1-delta)*kv(t-1);
            cv(t) = yv(t) - iv(t);
            rhseev_rbv(t) = exp( sum(x*solved_eta2_rhsee ) );
            rbv(t)    = 1/( cv(t) - psi*hv(t)^(1+1/nu)/(1+1/nu) ) / rhseev_rbv(t);
            rhseev(t) = rhseev_rbv(t) * rbv(t);
        end

        
        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next - erk_next) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;

    end
    % End of T period simulation

%% Pick up bank run events

window = 3;  % how many periods before and after runs to plot

% pick up all runV = 1

norun_num  = 0;
run_num    = 0;

sum_p_normal = 0;

rkhat_run  = zeros(1);
rkhats_run = zeros(1);
h_run   = zeros(1);
y_run   = zeros(1);
c_run   = zeros(1);
i_run   = zeros(1);
k_run   = zeros(1);
pib_run = zeros(1);
n_run   = zeros(1);
l_run   = zeros(1);
rb_run  = zeros(1);
pr_run  = zeros(1);
shock_run = zeros(1);
xi_run  = zeros(1);
tfp_run = zeros(1);

for t = t_drop:T-window

    if xiv(t) < 1

        run_num = run_num + 1;

        rkhat_run(run_num,1:1+2*window)  = rkhat(t-window:t+window);
        rkhats_run(run_num,1:1+2*window) = rkhats(t-window:t+window);
        h_run(run_num,1:1+2*window)   = hv(t-window:t+window);
        y_run(run_num,1:1+2*window)   = yv(t-window:t+window);
        c_run(run_num,1:1+2*window)   = cv(t-window:t+window);
        i_run(run_num,1:1+2*window)   = iv(t-window:t+window);
        k_run(run_num,1:1+2*window)   = kv(t-window:t+window);
        pib_run(run_num,1:1+2*window) = pibv(t-window:t+window);
        n_run(run_num,1:1+2*window)   = nv(t-window:t+window);
        l_run(run_num,1:1+2*window)   = lv(t-window:t+window);
        rb_run(run_num,1:1+2*window)  = rbv(t-window:t+window);
        pr_run(run_num,1:1+2*window)  = prv(t-window:t+window);
        shock_run(run_num,1:1+2*window) = shocks(t-window:t+window);
        tfp_run(run_num,1:1+2*window) = av(t-window:t+window);
        xi_run(run_num,1:1+2*window)  = xiv(t-window:t+window);

    else
        norun_num   = norun_num + 1;
    end

end

% BC moments
uncon_run_all = run_num / (T-window-t_drop)
avg_y_all  = mean(yv(t_drop:T));
avg_c_all  = mean(cv(t_drop:T));
avg_h_all  = mean(hv(t_drop:T));
avg_l_all  = mean(lv(t_drop:T));
avg_n_all  = mean(nv(t_drop:T));
avg_k_all  = mean(kv(t_drop:T));
avg_pr_all = mean(prv(t_drop:T));
avg_rb_all = mean(rbv(t_drop:T));

std_y_all  = std(yv(t_drop:T));
std_c_all  = std(cv(t_drop:T));
std_h_all  = std(hv(t_drop:T));
std_l_all  = std(lv(t_drop:T));
std_n_all  = std(nv(t_drop:T));
std_k_all  = std(kv(t_drop:T));
std_pr_all = std(prv(t_drop:T));
std_rb_all = std(rbv(t_drop:T));

% run dynamics
rkhat_run_avg  = mean(rkhat_run,1);
rkhats_run_avg = mean(rkhats_run,1);
h_run_avg   = mean(h_run,1);
y_run_avg   = mean(y_run,1);
c_run_avg   = mean(c_run,1);
i_run_avg   = mean(i_run,1);
k_run_avg   = mean(k_run,1);
pib_run_avg = mean(pib_run,1);
n_run_avg   = mean(n_run,1);
l_run_avg   = mean(l_run,1);
rb_run_avg  = mean(rb_run,1);
pr_run_avg  = mean(pr_run,1);
shock_run_avg = mean(shock_run,1);
xi_run_avg  = mean(xi_run,1);
tfp_run_avg  = mean(tfp_run,1);

%% pick up runV = 1 and no run some periods before

norun_period = 12;
run_num2 = 0;

rkhat_run2 = zeros(1);
rkhats_run2  = zeros(1);
h_run2   = zeros(1);
y_run2   = zeros(1);
c_run2   = zeros(1);
i_run2   = zeros(1);
k_run2   = zeros(1);
pib_run2 = zeros(1);
n_run2   = zeros(1);
l_run2   = zeros(1);
rb_run2  = zeros(1);
pr_run2  = zeros(1);
shock_run2 = zeros(1);
tfp_run2 = zeros(1);
xi_run2  = zeros(1);

norm_num = 0;
norm_nega = 0;
norm_y   = 0;
norm_c   = 0;
norm_l   = 0;
norm_k   = 0;
norm_n   = 0;
norm_pr  = 0;
norm_rb  = 0;
norm_h   = 0;
norm_pib = 0;
norm_posi_num = 0;
norm_l_posi = 0;

for t = t_drop:T-window

    if sum(runv(t-norun_period:t-1)) == 0
        
        if xiv(t) < 1

        run_num2 = run_num2 + 1;

        rkhat_run2(run_num2,1:1+2*window)  = rkhat(t-window:t+window);
        rkhats_run2(run_num2,1:1+2*window) = rkhats(t-window:t+window);
        h_run2(run_num2,1:1+2*window)   = hv(t-window:t+window);
        y_run2(run_num2,1:1+2*window)   = yv(t-window:t+window);
        c_run2(run_num2,1:1+2*window)   = cv(t-window:t+window);
        i_run2(run_num2,1:1+2*window)   = iv(t-window:t+window);
        k_run2(run_num2,1:1+2*window)   = kv(t-window:t+window);
        pib_run2(run_num2,1:1+2*window) = pibv(t-window:t+window);
        n_run2(run_num2,1:1+2*window)   = nv(t-window:t+window);
        l_run2(run_num2,1:1+2*window)   = lv(t-window:t+window);
        rb_run2(run_num2,1:1+2*window)  = rbv(t-window:t+window);
        pr_run2(run_num2,1:1+2*window)  = prv(t-window:t+window);
        shock_run2(run_num2,1:1+2*window) = shocks(t-window:t+window);
        tfp_run2(run_num2,1:1+2*window)  = av(t-window:t+window);
        xi_run2(run_num2,1:1+2*window)  = xiv(t-window:t+window);
        
        else
            
        norm_num = norm_num + 1;
        norm_y(norm_num)   = yv(t);
        norm_c(norm_num)   = cv(t);
        norm_l(norm_num)   = lv(t);
        norm_k(norm_num)   = kv(t);
        norm_n(norm_num)   = nv(t);
        norm_pr(norm_num)  = prv(t);
        norm_rb(norm_num)  = rbv(t);
        norm_h(norm_num)   = hv(t);
        norm_pib(norm_num) = pibv(t);

        if pibv(t) < 0
            norm_nega = norm_nega+1;
        else
            norm_posi_num = norm_posi_num + 1;
            norm_l_posi = norm_l_posi + lv(t);
        end
        
        end
    end
end

run_norun_before = run_num2 / (T-window-t_drop)

rkhat_run2_avg  = mean(rkhat_run2,1);
rkhats_run2_avg = mean(rkhats_run2,1);
h_run2_avg   = mean(h_run2,1);
y_run2_avg   = mean(y_run2,1);
c_run2_avg   = mean(c_run2,1);
i_run2_avg   = mean(i_run2,1);
k_run2_avg   = mean(k_run2,1);
pib_run2_avg = mean(pib_run2,1);
n_run2_avg   = mean(n_run2,1);
l_run2_avg   = mean(l_run2,1);
rb_run2_avg  = mean(rb_run2,1);
pr_run2_avg  = mean(pr_run2,1);
shock_run2_avg = mean(shock_run2,1);
tfp_run2_avg = mean(tfp_run2,1);
xi_run2_avg  = mean(xi_run2,1);

norm_y_avg   = mean(norm_y);
norm_c_avg   = mean(norm_c);
norm_l_avg   = mean(norm_l);
norm_k_avg   = mean(norm_k);
norm_n_avg   = mean(norm_n);
norm_pr_avg  = mean(norm_pr);
norm_rb_avg  = mean(norm_rb);
norm_h_avg   = mean(norm_h);
norm_pib_avg = mean(norm_pib);

nega_prob = norm_nega/norm_num;

norm_l_posi = norm_l_posi/norm_posi_num;

figure

subplot(2,3,1)
histogram(lv(t_drop:T),'Normalization','probability')
hold on
plot([l_stss_nopol,l_stss_nopol],ylim,'r--','Linewidth',3);
hold on
plot([thre_l,thre_l],ylim,'k:','Linewidth',3);
set(gca,'FontSize',20);
grid on
title({'leverage','all time'})

subplot(2,3,2)
histogram(nv(t_drop:T),'Normalization','probability')
hold on
plot([n_stss_nopol,n_stss_nopol],ylim,'r--','Linewidth',3);
%hold on
%plot([thre_n,thre_n],ylim,'k:','Linewidth',3);
set(gca,'FontSize',20);
title({'net worth','all time'})

subplot(2,3,3)
histogram(rbv(t_drop:T),'Normalization','probability')
hold on
plot([rb_stss_nopol,rb_stss_nopol],ylim,'r--','Linewidth',3);
set(gca,'FontSize',20);
title({'deposit rate','all time'})

subplot(2,3,4)
histogram(norm_l,'Normalization','probability')
hold on
plot([l_stss_nopol,l_stss_nopol],ylim,'r--','Linewidth',3);
hold on
plot([thre_l,thre_l],ylim,'k:','Linewidth',3);
set(gca,'FontSize',20);
grid on
title({'leverage',['no run for ',num2str(norun_period),' periods']})

subplot(2,3,5)
histogram(norm_n,'Normalization','probability')
hold on
plot([n_stss_nopol,n_stss_nopol],ylim,'r--','Linewidth',3);
%hold on
%plot([thre_n,thre_n],ylim,'k:','Linewidth',3);
set(gca,'FontSize',20);
title({'net worth',['no run for ',num2str(norun_period),' periods']})

subplot(2,3,6)
histogram(norm_rb,'Normalization','probability')
hold on
plot([rb_stss_nopol,rb_stss_nopol],ylim,'r--','Linewidth',3);
set(gca,'FontSize',20);
title({'deposit rate',['no run for ',num2str(norun_period),' periods']})
